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In this paper, I present a precise Quantum Monte Carlo calculation at 
finite temperature for a very large number (many thousands) of bosons in a 
harmonic trap, which may be anisotropic. The calculation applies directly 
to the recent experiments of Bose-Einstein condensation of atomic vapors in 
magnetic traps. I show that the critical temperature of the system decreases 
with the interaction. I also present profiles for the overall density and the one 
of condensed particles, and obtain excellent agreement with solutions of the 
Gross-Pitaevskii equation. 
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The achievement of Bose-Einstein condensation (BEC) in dilute atomic vapors of 87 Rb 
TJ and 23 Na J|] has created extraordinary experimental and theoretical activity. The exper- 
iments pose fundamental theoretical questions, ranging from a clear understanding of the 
dynamics of BEC in a finite sample || to the behavior under time-dependent perturbations. 

Besides these tantalizing and very controversial questions about the dynamics of the 
quantum coherence, there are many questions about the statics (thermodynamics) of the 
BE-condensate in a harmonic trap. In this paper, I show that we can obtain the static 
properties from Path-Integral Quantum Monte Carlo (QMC) calculations which suffer from 
no systematic uncertainties other than the choice of the interatomic (pseudo-)potential. I 
will present results for 10000 particles in isotropic and anisotropic traps around the BEC 
transition point, i.e. at the temperatures and particle number of current experimental 
interest. 

The hamiltonian of the system is given by 

H = f> 2 /2m + £ V(nj) + f (coy + coy + coy) (1) 

i=l i,j 

where V is the interatomic potential M. This potential has many bound states, which, in 
the experiment, are inaccessible due to the kinematics of two-particle collisions. In ther- 
modynamic calculations, as the one presented here, it is necessary to work with a pseudo- 
potential. I simply adopt a hard-core term of radius - as others have done before. More 
refined choices can easily be implemented. 

The physical parameters - and the units - are as in ref. ||: we consider two harmonic 
traps (isotropic and anisotropic), and units in which all frequencies are equal to 1 in the 
isotropic case, and co x = co y = 1, co z = \/8 in the anisotropic case. The anisotropy seems to 
be close to the experimental situation. The value a for the hard-core radius is taken to be 
ao = 0.0043 in the above units 0. 

To situate the computations of the interacting case, it is useful to remember the basic 
result for non-interacting bosons (in our isotropic potential, and with the above units) || 
0: for large N, the transition temperature T c scales as T c ~ T C N 1 ^ with T c ~ 0.94. For 
iV = 10.000 particles this gives (3 C ~ 0.05. Let us mention in passing that at this value of 
N there are no more detectable differences between the canonical and the grand-canonical 
formalism [§]. 

The partition function Z, the trace over the density matrix, satisfies the usual convolution 
equation: 

Z = J^ Hj dR p(R, R P ,P) = ^ E / • • • / dRdR 2 ■ ■ ■ dR M p(R, R 2 , r) . . . p(R M , R P , r) 

(2) 

Here r = /3/M; R is the 3iV-dimensional vector of the particle coordinates R = 
(r 1; r 2 , . . . ,r N ), and R p denotes the vector with permuted particle labels: R p = 
(r P (i),r P (2), ■ ■ ■ ,r P ( N )). 

In the limit M — > oo, the exact convolution equation eq. (0) reduces to the Feynman- 
Kac formula (cf [|IIJ])> as a consequence of the Trotter break-up exp[r(Hi + H 2 )] ~ 
exp[riJi] expfri^] for r — ► 0. 
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However, it is the basic technical lesson of past QMC work [|TT| |~2] that eq. (0) should 
only be taken to intermediate temperatures r such that the approximation of the complete 
density matrix p(R, R', r) by a product of pair terms becomes accurate 



p(R, R', r) ~ n Pifc, r>, r) l[ P^^l , (3) 

The one-particle density matrices pi(rj,r 4 ',r) in eq. (^|) describe a particle in the harmonic 
potential: an exact solution is known (cf., e.g. JTOf); and can easily be sampled. This 
product of one-particle density matrices is used as a priori probability density in the sense 
of ref. [0. A very efficient method is used to generate a sequence of configurations, and 
to sample permutations. The second part in eq. (|3|) is equal to 1 for V = 0. It plays the 
role of a correction term due to two-body interaction. The Metropolis algorithm is used 
to incorporate this term, and to generate configurations which are distributed according to 
eq. (D. 

To obtain the correction factor we need the two-particle density matrix. Fortunately, for 
the special case of the harmonic potential, the two-particle hamiltonian - and therefore the 
density matrix - decouples into a center-of-mass term (which is unaffected by the interaction) 
and a term describing the relative motion. The latter is given by: 

H rel = p*/2p + V(r) + + uy + ^V) (4) 

where p = m/2 is the reduced mass. 

The eigenfunctions of the pure hard-sphere potential (eq. without the terms in u>) can 
be obtained exactly ( |fL3|). Adding an isotropic trap potential, the same may still be possible. 
However, since I am interested in treating the general anisotropic trap, I have obtained the 
relative-motion density matrix from the exact solution of the hard-core potential density 
matrix ph c , into which I incorporate the harmonic term via the Trotter break-up. 

p2,rd{r, r', t) ~ X(r) x p hc (r, r') x X{r') (5) 

with X{r) = exp[— rpiuj^x 1 + uj^y 1 + uj 2 z z 2 )/4\. The r finally retained in the simulation 
has to be chosen to accommodate both eq. (0) and eq. (H). Extensive tests have convinced 
me that for the density attained, a value of r = 0.01 is appropriate. This value is about 
6 orders of magnitude larger than the one which would have to be used with the simple 
Trotter break-up. 

For the value of r = 0.01 which was found sufficient, we can explicitly compute the effec- 
tive range of the interaction, beyond which the correction factor in eq. @ is practically one. 
This range turns out to be about 0.2. It is evident that one may introduce a 3— dimensional 
grid, with the grid size larger than the interaction range. At any time, particles are assigned 
to boxes formed by the grid. In evaluating the pair density matrix eq. (0), an efficient al- 
gorithm can be set up to compute the correction term only for close-by pairs. As a result, 
the program behaves gracefully as is increased, and the actual limit of the calculation is 
rather given by the memory requirements than by CPU considerations. 

The program has been checked very carefully. In the non-interacting case, I am able 
to reproduce all the exactly known results. For an isotropic trap, and in the absence of 
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interactions, e. g. I have reproduced the analytically known results for the condensate 
fraction N /N, which is plotted in fig. 1. I also plot the same quantity with interactions: 
the condensate fraction clearly decreases with respect to the noninteracting case. From this 
I conclude that T c decreases. Notice that the interaction influences the condensate fraction 
much more than the finite-size effects for the non-interacting gas, which are also shown. 

On the left side of fig. 2, I show the corresponding density profiles p(x) = 
J J dydzp(x, y, z) of the particles. The center density increases sharply as T is lowered. 
This is the hallmark of real-space BEC in the confining potential. 

We pause for a moment to discuss BEC in the path-integral framework: below T c , par- 
ticles have a finite probability to belong to extended permutation cycles of length I. Notice 
that in the confined geometry the off-diagonal one-particle density matrix at large separa- 
tions trivially vanishes, instead of going to a value proportional to Nq (as in the translational 
invariant system). In the present simulation, I instead obtain the number of condensed par- 
ticles from the permutations of the system: the maximum length of / which has non-zero 
probability is equal to iVo- For non-interacting particles in the limit T — > 0, the probability 
to belong to extended permutation cycles is independent of I (having all particles in the 
condensate does therefore not mean that they are all on the same cycle). At any temper- 
ature, noninteracting bosons which are on the different permutation cycles are statistically 
independent |12|] - a property which is used in the QMC code to generate the a priori prob- 



ability. The spatial distribution of particles in a cycle of length I is given by the diagonal 
density matrix at temperature 1(3 [ |10|] : 



p(x, x, 1(3) ~ exp(— mux 2 tanh — (31) (6) 

Since (3 C ~ 1/iV 1 / 3 , the particles on long cycles with I >> iV 1//3 (especially / ~ N) are 
distributed according to the lowest single-particle state ^l(x) ~ exp(— mux 2 ). 

As we introduce interactions, we have to give up the concept of condensation into the 
lowest single-particle state. However, we preserve the two other essential features: below the 
transition, long permutation appear, and particles on long permutation cycles are distributed 
identically: they populate the macroscopic quantum state. We can thus gain access at 
the distribution of the condensed particles by computing the density profiles as before, but 
restricted to particles on cycles longer than some value l m in- I have verified explicitly that for 
l min ~ 20 the density profile does not depend on this parameter, i. e. that the condensation 
concerns the ground state. The density profiles for the condensed particles is given on the 
right side of fig. (2); the distributions are normalized to 1 at any temperature. As we 
expect, the distribution of condensed particles is broader at low temperature, where the 
many particles in the condensate strongly repel each other. 

Can we obtain a more quantitative description of the condensate at finite temperature? 
To answer the question, I have computed the solution of the isotropic Gross-Pitaevskii (G-P) 
equation [|14[] \TE\ for the same value of a® as is used in the QMC-calculation, and for the 



number of particles corresponding approximately to the condensate fraction, as obtained 
from fig. 1. The results ||16| are plotted together with p(x) for the condensed particles on 
the right side of fig. 2. The agreement is truly remarkable for all the three curves, especially 
since there are no adjustable parameters. Even the number of condensed particles has been 
simply taken from fig. 1, without trying to optimize the fit. It can thus be said that, to a 
very high precision, the condensed particles are described by the G-P wavefunction. 
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The numerical results presented in fig. 2 suggest the following quantitatively correct 
picture for the Bose-Einstein condensation in a trap at finite temperature: below the 
(interaction-dependent) T c , N (T) particles are condensed. To a very high precision, these 
particles are distributed according to the appropriate G-P wavefunction. The non-condensed 
particles are distributed as in the free case and are very much spread out. There is very 
little mutual interaction between condensed and non-condensed particles: on the one hand, 
the non-condensed particles are very dilute in the central region of the trap, where most of 
the particles present are "condensed" , on the other hand the G-P wavefunction disturbs the 
non-condensed particles only on a small portion of their support. 

The above picture allows us to understand the competition of energy gain end entropy de- 
struction which underlies the condensation into the macroscopic wavefunction. The balance 
of entropy is the same for the non-interacting and interacting case, since the non-condensed 
part is undisturbed by the interaction, and the condensed one has zero entropy. The en- 
ergy of the G-P wavefunction is of course much higher than in the absence of interactions. 
Condensation is therefore less favored, and we understand that the critical temperature 
must decrease with the interaction, as shown in fig. 1. It will be very simple to perform 
a one-parameter (in No) variational minimization of the free energy which describes the 
non-condensed particles as free, and which uses the energy of the G-P wavefunction. 

Rather than to proceed in this direction, I finally present some results for the anisotropic 
trap, which is the case of direct interest to experiments. In fig. 3, I present density plots 
for this case at different temperatures, and compare to the solution of the anisotropic G-P 
wavefunction with N = 10000 ||. As before, there is very close agreement between the 
two calculations. The calculation also confirms that, in the limit of zero temperature, the 
number of condensed particles, N , seems to be extremely close to N. This agrees with 
recent work in which the interaction of the G-P groundstate with non-condensed particles 



was studied in the Bogolubov approximation [17]. I have performed the same calculation as 



in fig. 1 for the isotropic case (the noninteracting curve is easily computed, even for finite 
iV H), and the results are analogous. For example, at /3 — 0.06, the noninteracting gas in 
the anisotropic trap has N /N = 0.76, but I only find a value of 0.6 for the interacting case. 

Even though the main impetus in this paper has been on the dependence of the critical 
temperature on the interaction, and on the comparison of the condensate wave function with 
the Gross-Pitaevskii formalism, it should be evident that the QMC calculation goes much 
beyond these results. We can obtain complete information on inter-particle correlation, 
compute the complete thermodynamics of the system, etc. In particular, it is possible |T8 



to obtain the nonclassical moment of inertia which has attracted some attention lately |19 



and which will be of experimental relevance as soon as rotating traps will become available. 
If the above picture is correct, this moment (which basically follows the normal part of the 
gas) should depend on the interaction only through Nq. 

I would also like to draw the reader's attention to the fact that the very efficient algo- 
rithm is of general usefulness for the study of weakly interacting bosons. Other possible 
applications are bosons in 3 dimensions without confining potential, but particularly in 2 
dimensions, where the effect of a small interaction will be even more important than in the 
case treated here, since the non-interacting gas has no phase transition at all. To stimulate 
work in this well as to facilitate direct comparison with experiences on trapped 

bosons, I will make available the FORTRAN code used in the present investigation. 
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Figure Captions 

1. Ratio of condensed particles N /N vs reduced temperature T = TN 1 ^ 3 in an isotropic 
trap with ou — 1 for the noninteracting case in the thermodynamic limit (N — > oo, full 
line) and for N = 10000 (dashed line), as well as for 10000 interacting particles with 
a = 0.0043. The number of condensed particles decreases with the interaction. 

2. Density profile vs x for temperatures (3 = 0.06(f = 0.78), f3 = 0.07(T = 0.66), 
(3 = 0.12(T = 0.39) (left side, from below). On the right side is plotted the density 
profile of the "condensed" particles for the same temperatures (from above, the curves 
are normalized to one). These density profiles are in excellent agreement with the 
Gross-Pitaevskii solutions for the iVo obtained from fig. 1 (dotted lines). The values 
used are N = 2000 for (3 = 0.06 (upper) JV = 4000 for (3 = 0.07 (middle) iVo = 8000 
for = 0.12 (lower). 

3. Density profile vs x for an anisotropic trap with iV = 10000. As in fig. 2, we plot 
the total density on the left, and the normalized density of condensed particles on the 
right. The temperatures are (3 = 0.06 (dotted line), and (3 = 0.16 (full line). The 
dashed curve is the Gross-Pitaevskii solution for iV = 10000. 
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